home *** CD-ROM | disk | FTP | other *** search
/ Turnbull China Bikeride / Turnbull China Bikeride - Disc 1.iso / ARGONET / PD / MATHS / RLAB / RLAB125.ZIP / !RLaB / toolbox / roots < prev    next >
Text File  |  1995-03-05  |  4KB  |  166 lines
  1. //-------------------------------------------------------------------//
  2. //  roots
  3.  
  4. //  Syntax: roots(P)
  5.  
  6. //  Description:
  7.  
  8. //  Find roots of polynomial P of degree n.
  9. //
  10. //  Case 1:
  11. //    P is a vector (row or column) with n+1 components representing
  12. //    polynomial coefficients in descending powers.
  13. // 
  14. //  Case 2:
  15. //    P is a linked-list of n+1 square matrices, the n+1 matrices in P 
  16. //    representing polynomial coefficients in descending powers.
  17. //
  18. //  Return a column vector containing roots of P.
  19.  
  20. //  Example 1:
  21. //  The following script finds the roots of 
  22. //  x^3 + 3*x^2 + 5 = 0
  23. // 
  24. //  > p=[1,3,0,5]
  25. //   p =
  26. //          1          3          0          5
  27. //  > r=roots(p)
  28. //   r =
  29. //             0.213 - 1.19i
  30. //             0.213 + 1.19i
  31. //                -3.43 + 0i
  32. //  > q=<<1;3;0;5>>;
  33. //  > r=roots(q)
  34. //   r =
  35. //             0.213 - 1.19i
  36. //             0.213 + 1.19i
  37. //                -3.43 + 0i
  38. //  
  39. //  Example 2:
  40. //  Find the damped and undamped frequencies of the following oscillation 
  41. //  system:
  42. //           2
  43. //          d u     du
  44. //        M --- + C -- + K u = 0
  45. //            2     dt
  46. //          dt  
  47. //
  48. //  where M =   5  0.5      C =  5   1      K = 50   10
  49. //            0.5    5           1   5          10   50
  50. //
  51. //  The characteristic equation corresponding to the governing equation
  52. //  is
  53. //        M*s^2 + C*s + K = 0
  54. //
  55. //  > M=[5,0.5;0.5,5]; C=[5,1;1,5]; K=[50,10;10,50];
  56. //  > // undamped frequencies
  57. //  > uf = roots(<<M;zeros(2,2);K>>)
  58. //  >  uf =
  59. //  >        -3.97e-16 - 2.98i
  60. //  >        -3.97e-16 + 2.98i
  61. //  >          2.84e-16 - 3.3i
  62. //  >          2.84e-16 + 3.3i
  63. //  The undamped vibration frequencies are 2.98 and 3.3 rad/s.
  64. //  > // damped frequencies
  65. //  > df = roots(<<M;C;K>>)
  66. //  >  df =
  67. //  >           -0.444 - 2.95i
  68. //  >           -0.444 + 2.95i
  69. //  >           -0.545 - 3.26i
  70. //  >           -0.545 + 3.26i
  71. //  So we know the damped vibration frequencies are 2.95 and 3.26 rad/s.
  72. //
  73. //  Example 3:
  74. //  Same as Example 2, but  C = 35   5
  75. //                               5  35
  76. //                            
  77. //  > M=[5,0.5;0.5,5]; C=[35,5;5,35]; K=[50,10;10,50];
  78. //  > f = roots(<<M;C;K>>)
  79. //  >  f =
  80. //        -5.16
  81. //        -4.82
  82. //        -2.12
  83. //        -1.84
  84. //  All roots are negative real ===> no oscillation but stable. 
  85.  
  86. //  See Also: poly
  87.  
  88. //  Tzong-Shuoh Yang 
  89. // (tsyang@ce.berkeley.edu)  5/7/94  scalar polynomial roots
  90. //                           3/2/95  matrix polynomial roots
  91. //
  92. //-------------------------------------------------------------------//
  93. roots = function(P)
  94. {
  95.    local(P);
  96.  
  97.    if (class(P) == "num") {
  98.      // scalar polynomial
  99.      if (P.n <= 1) {
  100.         error("Argument of roots() is a scalar, not a polynomial.");
  101.      else if (min(size(P)) != 1) {
  102.         error("Argument of roots() must be a column or a row vector.");
  103.      }}
  104.      p = P[:]';
  105.      // trim leading and trailing zeros
  106.      z  = find(p != 0);
  107.      nz = p.n - z[z.n];
  108.      p  = p[z[1]:z[z.n]];
  109.      // find companion matrix
  110.      c = compan(p);
  111.    else if (class(P) == "list") {
  112.      // matrix polynomial
  113.      n = max(strtod(members(P)));
  114.      if (n <= 1) {
  115.         error("Argument list of roots() is not a matrix polynomial.");
  116.      }
  117.      ns = P.[1].nr;
  118.      if (ns != P.[1].nc) {
  119.         error("Argument list components of roots() must be square matrices.");
  120.      }
  121.      for (i in 2:n) {
  122.         if (P.[i].nr != ns || P.[i].nc != ns) {
  123.         error("Argument list components of roots() must be the same size.");        
  124.         }
  125.      }
  126.      // find leading and trailing zero matrices
  127.      z = [];
  128.      for ( i in 1:n) {
  129.        if (any((P.[i] != 0)[:])) {
  130.           z = [z,i];
  131.        }
  132.      }
  133.      nz = (n - z[z.n])*ns;     
  134.      // find companion matrix c
  135.      n0 = z[1];
  136.      n1 = z[z.n];
  137.      P.[n0] = factor(P.[n0],"g");
  138.      for (i in (n0+1):n1) {
  139.          P.[i] = -backsub(P.[n0], P.[i]);
  140.      }
  141.      n2 = (n1-n0)*ns;
  142.      c = zeros(n2,n2);
  143.      for (i in 1:(n2-ns)) {
  144.         c[i;ns+i] = 1;
  145.      }
  146.      j = (n2-ns+1):n2;
  147.      k = j;
  148.      for(i in (n0+1):n1) {
  149.         c[j;k] = P.[i];
  150.         k = k - ns;
  151.      }
  152.    else {
  153.      error("Wrong argument class in roots().");
  154.    }}}   
  155.    // find roots and stuff with trailing zero(s)
  156.    r = [zeros(nz,1);eign(c).val'];   
  157.    // trim complex part if all real
  158.    if (all(imag(r) == 0)) {
  159.       r = real(r);
  160.    }
  161.    // if you don't want the roots sorted in ascending order, 
  162.    // comment out the next line
  163.    r = sort(r).val;
  164.    return r;
  165. };
  166.